↜ Back to index Introduction to Numerical Analysis 2
Lecture a8
Survey
!!! Please fill out the class survey 授業改善のための学生アンケート.
アカンサスポータル＞学務情報サービス＞アンケート回答＞授業アンケート
LU decomposition
The Gaussian elimination algorithm that we implemented proceeds in two steps:
Convert the system Ax = b into a system Ux = \tilde b with a upper triangular matrix U and a modified vector \tilde b. This performs elimination by adding rows of the matrix and swapping the rows for pivoting.
Solve the system Ux = \tilde b for x.
Let us try to understand how to compute \tilde b from b by keeping track of the changes to b during the run of the Gaussian elimination algorithm.
No pivoting
To simplify, we do not consider pivoting first.
For column 1, the algorithm performs the following:
a_{ij} \leftarrow a_{ij} - \frac{a_{i1}}{a_{11}} a_{1j}, i = 2, \ldots, n, j = 1, ..., n.
b_i \leftarrow b_i - \frac{a_{i1}}{a_{11}} b_1, i = 2, \ldots, n.
Let us set l_{i1} := \frac{a_{i1}}{a_{11}}, i > 1.
Similarly for column 2:
a_{ij} \leftarrow a_{ij} - \frac{a_{i2}}{a_{22}} a_{2j}, i = 3, \ldots, n, j = 2, ..., n.
b_i \leftarrow b_i - \frac{a_{i2}}{a_{22}} b_2, i = 3, \ldots, n.
Let us again set l_{i2} := \frac{a_{i2}}{a_{22}}, i > 2. Note that a_{i2} and a_{22} are not the original components of A, but the new values that were assigned to them when working on column 1.
The algorithm continues like this for columns i = 3, \ldots n.
We therefore see that
b_1 does not change and so \tilde b_1 = b_1.
b_2 changes only when working on column 1, and the final value is
\tilde b_2 = b_2 - l_{21} \tilde b_1
b_3 changes only when working on columns 1 and 2, first b_3 \leftarrow b_3 - l_{31} b_1 and then b_3 \leftarrow b_3 - l_{32} b_2, where b_2 is the value that we assigned when working on column 1, that is, it is \tilde b_2.
Therefore the final value of b_3 happens to be
\tilde b_3 = b_3 - l_{31} \tilde b_1 - l_{32} \tilde b_2.
We can continue like this, recovering the pattern \tilde b_i = b_i - l_{i1} \tilde b_1 - \cdots - l_{i,i-1} \tilde b_{i-1}.
We saw this pattern before in Lecture 4, when we were solving a linear system with a lower triangular matrix L = \begin{pmatrix} 1 & 0 &&\cdots & 0\\ l_{21} & 1 & 0 && \vdots\\ l_{31} & l_{32} & 1 & 0\\ \vdots& & & \ddots &0\\ l_{n1} & l_{n2} & \cdots & l_{n,{n-1}} & 1\\ \end{pmatrix}
We can therefore find \tilde b by solving the linear system L \tilde b = b.
All we need to do is to store the values l_{ij} obtained during the Gaussian elimination in a matrix L.
The resulting two matrices L and U are an LU decomposition of the matrix A. It is not difficult to see that A = L U (the matrix product).
Pivoting
Unfortunately, not every regular matrix has an LU decomposition. Sometimes we need to swap the rows during the Gaussian elimination algorithm.
When considering column k of A, we might need to swap the row k with some row m, k < m \leq n, and also swap b_k and b_m. To solve U x= \tilde b correctly, the meaning of \tilde b_k and \tilde b_m must be also swapped.
Consider the matrix A = \begin{pmatrix} 0 &1\\ 1 & 0 \end{pmatrix}
Clearly, when processing the first column we need to swap the rows since a_{11} is 0, and therefore the result of the algorithm becomes the combined LU matrix \begin{pmatrix} 1 &0\\ 0 & 1 \end{pmatrix} which represents L = \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix} \qquad U = \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix}.
If we use it to solve A x = b, we get L \tilde b = b \quad \Rightarrow \quad \tilde b = b and U x = \tilde b \quad \Rightarrow \quad x = b
But that is wrong since the actual solution of Ax = b is x = (b_2, b_1).
It turns out that keep track of the swaps is very simple.
We need to swap the rows k and m of L at the same time we are swapping the rows of A.
To keep track of the swaps of the components of b, we use the following idea:
We initialize an array p_1, …, p_n with values p_i = i, so that b_{p_i} = b_i initially.
But then every time we swap the rows k and m of A and L, we also swap the entries p_k and p_m.
Then when we want to find \tilde b from b, we solve L \tilde b = b_p, where b_p is the vector with components b_{p_1}, …, b_{p_n}.
This gives an LUP decomposition.
Note that it is OK for the matrix U to have zero values on the main diagonal.
Consider the matrix from the previous example: A = \begin{pmatrix} 0 &1\\ 1 & 0 \end{pmatrix}
We start with array p with values (1,2). When the rows are swapped, we also swap the values p_1 and p_2 and therefore we end up with an LUP decomposition
L = \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix} \qquad U = \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix} \qquad p = (2, 1).
If we use it to solve A x = b, we start with reordering the right hand side vector b according to the values in p: b_p := (b_{p_1}, b_{p_2}) = (b_2, b_1), and then we get L \tilde b = b_p \quad \Rightarrow \quad \tilde b = b_p and U x = \tilde b \quad \Rightarrow \quad x = b_p
This is the correct result x = (b_2, b_1)!
Storing L
Now that we know how to construct U, L and p, let us think about storing them in memory.
The array p is simply a one-dimensional array of length n.
As for matrix L, we could store it in a two-dimensional array, separate for the array a we use for A and U. However, there is a better way.
Notice that all the components of U below the diagonal are always 0, and therefore we do not really need to store their value. At the same time, all the interesting components of L are under the diagonal, because the only possible values on the diagonal of L are 1. We therefore might store both U and L in the same array a. Furthermore, when we perform swapping of rows, we need to do them for both matrices at the same time so this way of storing them makes it more convenient. We end up with an array that represents the matrix
\begin{pmatrix} u_{11} & u_{12} &&\cdots & u_{1n}\\ l_{21} & u_{22} & u_{23} && \vdots\\ l_{31} & l_{32} & u_{33} & u_{34}\\ \vdots& & & \ddots &u_{n-1,n}\\ l_{n1} & l_{n2} & \cdots & l_{n,{n-1}} & u_{nn}\\ \end{pmatrix} = L + U - I, where I is the identity matrix.
Example
Suppose that the input array contains \begin{pmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{pmatrix}
After processing the first column, if no pivoting is necessary, we obtain the array \begin{pmatrix} a_{11} & a_{12} & a_{13}\\ \frac{a_{21}}{a_{11}} & a_{22} - \frac{a_{21}}{a_{11}} a_{12} & a_{23} - \frac{a_{21}}{a_{11}} a_{13} \\ \frac{a_{31}}{a_{11}} & a_{32}- \frac{a_{31}}{a_{11}}a_{12} & a_{33}- \frac{a_{31}}{a_{11}} a_{13} \end{pmatrix}
The algorithm
The algorithm for LU decomposition with pivoting can be written as follows.
Initialize the two dimensional array a by reading the elements of matrix A.
Initialize an integer array p with values p_1 = 1, …, p_n = n.
For each column k = 1, \ldots, n-1 do
Find a pivot: find m satisfying k \leq m \leq n so that |a_{mk}| is the largest.
If |a_{m k}| < \varepsilon, continue with the next k.
If m \neq k, swap the row k and m of a, and swap p_k and p_m
Store a_{ik} \leftarrow \frac{a_{ik}}{a_{kk}}, \qquad i = k+1, \dots, n
Add rows (using the new values a_{ik}!) a_{ij} \leftarrow a_{ij} - a_{ik} a_{kj}, \qquad i = k+1, \ldots, n, \ j = k+1, \ldots, n Be careful about the range of j.
Repeat for next k.
The parameter \varepsilon is a small positive constant, for example, \varepsilon = 0.00001.
Using the LUP decomposition to solve the linear system
Once we have the LUP decomposition L, U and p for a matrix A, it is easy to quickly1 solve the linear system Ax = b.
The algorithm is:
Reorder the components of b using p and produce b_p. b_p has components (b_{p_1}, b_{p_2}, \ldots, b_{p_n}), where b_{p_1} is the component of b with index p_1, etc.
Solve L \tilde b = b_p for \tilde b.
Solve U x = \tilde b for x.
When should you use the LUP decomposition?
If you need to solve Ax = b for a given matrix A once or twice, just use the Gaussian elimination.
However, in practice, you might need to solve the system Ax = b with the same matrix A for many different right hand sides b. In that case, first compute the LUP decomposition and then use it to solve Ax = b much faster by solving LU x = b_p for each of the right hand sides b.
Note that for some matrices there are faster methods to find the solution. We will talk about one of them, the conjugate gradient method next time.
Report 2
Upload the code for the short report to Acanthus portal by Dec 3 (Fri). Maximum score is 150 points.
Make sure that the program compiles before submission. If the code does not compile, the score is automatically 0.
Make sure that there are no “index out of bounds” errors by compiling with
-fcheck=boundsoption as explained in Lecture 3.Make sure that the program produces the correct output for the example inputs. Due to numerical errors the values might be a bit different, but they should still be close to the ones given in the examples below.
More or fewer spaces or new lines are OK, but there must be no extra text in the output!
Do not submit the Emacs temporary file ending with
~!
Problem 1. Implement the LUP decomposition algorithm with pivoting explained above in Fortran.
Input: Number n followed by n^2 components of the matrix A in the same order as in the previous problems.
Output: The result of the LUP decomposition in the form of the combined LU matrix (the nonzero components of U are on the diagonal and above, the nontrivial components of L are below the diagonal) followed by the elements of the array p.
Download the test files (they are the same as in Lecture 6):
Expected outputs with the individual test files:
$ ./a.exe < data.txt
3.00000000 11.0000000 6.00000000
0.333333343 -2.66666675 -3.00000000
0.333333343 0.250000030 -0.249999881
3
2
1
$ ./a.exe < data2.txt
4.00000000 5.00000000 6.00000000 5.00000000 3.00000000
0.250000000 3.75000000 2.50000000 4.75000000 1.25000000
0.750000000 0.866666675 -1.66666675 -3.86666679 0.666666627
1.00000000 -0.533333361 -0.800000012 57.4399986 0.199999988
0.500000000 0.933333337 0.199999943 -8.98328796E-02 10.2179661
4
1
3
5
2
$ ./a.exe < data3.txt
1.00000000 0.00000000
0.00000000 1.00000000
2
1
$ ./a.exe < data4.txt
3.00000000 11.0000000 5.00000000
0.333333343 -2.66666675 -2.66666675
0.333333343 0.250000030 0.00000000
3
2
1Problem 2. Write a program that uses an LUP decomposition to solve A x = b.
Input: Number n, followed by n^2 components of the combined matrix LU, then n components of the array p and finally n components of the vector b.
Output: n components of the solution x_1, …, x_n, or an error message if U has a zero on the diagonal.
Download the new test files:
Expected outputs with the individual test files:
$ ./a.exe < sys1.txt
-4.99999952
4.99999952
-0.999999523
$ ./a.exe < sys2.txt
Error: matrix U has zero on the main diagonal
$ ./a.exe < sys3.txt
1.63013709
-0.188509524
3.41443419E-02
1.16540594E-02
5.29544093E-02Problem 3. Write a program that given the LUP decomposition prints out the matrix A.
First multiply L and U as matrices, and then reorder the rows of the product using p.
Input: Number n, followed by n^2 components of the combined matrix LU, then n components of the array p.
Output: n^2 components of the matrix A.
Use the same files as in Problem 2.
$ ./a.exe < sys1.txt
1.00000000 3.00000000 1.00000000
1.00000000 1.00000000 -1.00000000
3.00000000 11.0000000 6.00000000
$ ./a.exe < sys2.txt
1.00000000 3.00000000 1.00000000
1.00000000 1.00000000 -1.00000000
3.00000000 11.0000000 5.00000000
$ ./a.exe < sys3.txt
1.00000000 5.00000000 4.00000000 6.00000000 2.00000000
2.00000000 6.00000000 5.00000000 1.00000000 13.0000000
3.00000000 7.00000000 5.00000000 4.00000000 4.00000000
4.00000000 5.00000000 6.00000000 5.00000000 3.00000000
4.00000000 3.00000000 6.00000000 63.0000000 2.00000000 Hints.
The matrix product C = LU has the components:
c_{ij} = \sum_{k=1}^n l_{ik} u_{kj}
Using the fact that l_{ii} = 1 for any i and l_{ij} = 0 for j>i, and that u_{ij} = 0 for j < i, we get the formula c_{ij} = \begin{cases} u_{ij} + \sum_{k=1}^{i-1} l_{ik} u_{kj}, & i \leq j,\\ \sum_{k=1}^j l_{ik} u_{kj}, & i > j. \end{cases}
The row permutation is the “inverse” of the permutation of b we were using before:
a_{p_i, j} = c_{i, j}
in a time proportional to n^2↩︎